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INTRODUCTION 

Recent impetus, for research in unsteady separated flows stems from a wide range 

of applications from low- to high- Reynolds , number , Re. The physics of high-Re 

flows, in general, is quite complex and often involves multiple nonuniqueness and 

chaos, beyond simple unsteady separation. For the low-Re case, e.g. in the 

manuevering of fighter aircraft at high angle-of-attack in near- and post-stall 

regime, the vortex interaction dominates the flow field. The passage of vortices 

over the suction surface and their subsequent shedding leads to self-excited 

persistently unsteady flows. This flow field is extremely complicated due to the 

global effect of unsteady separated flow, coupled with the presence of hydrodynamic 

instabilities which may trigger transition and eventually lead to chaos. Besides 

supermaneuverability, interest also lies in this low-Re case because of the need for 
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design of efficient airfoil sections for Re in the range of 10-10 , for improving 

the performance of mini-RPV's (remotely piloted vehicles) operating at low 

altitudes, jet engine compressor and turbine blades, helicopter rotor blades, etc. 

For low-speed viscous flow without body forces, Re is the key similarity 

parameter. For flow over lifting bodies, the flow pattern is unique and steady for 

Re < Re , where Re r Is the value of Re at which transition first occurs. Near and 
cr’ . 

beyond Re cr , the flow is highly unsteady and it becomes imperative that this 
unsteady flow be better understood. Recent developments in the dynamical theory of 
low-dimensional nonlinear systems have provided a new and stimulating viewpoint 
concerning the onset of turbulence, as computer simulations of three or more coupled 
nonlinear first-order ordinary differential equations have led to chaotic solutions. 
Such an observation led Ruelle and Takens (1971) to hypothesize that transition to 
turbulence can be quantitatively explained by deterministic equations. Their work 
has produced exciting and profound results and the resulting theory has become 
widely known as the Ruelle-Takens theory of turbulence. They showed mathematically 
that, for a nonlinear initial-boundary value problem, with a large critical 
parameter such as Re, chaos resulted from repeated bifurcations of the solution. 
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The chaotic solutions were referred to as "Strange Attractors", to distinguish them 
from other ordinary attractors such as fixed points and limit cycles, i.e. steady- 
state and periodic solutions, respectively. The strange attractors are geometric 
entities in the phase or state space of the governing differential equations. In 
this context, fixed points are zero-dimensional entities, limit cycles are one- 
dimensional entities in the sense that they can be parameterized by a single degree 
of freedom, and strange attractors are more-than-one dimensional, perhaps even of 
fractional dimension. The important point is that the dimension of a strange 
attractor is often quite low. Indeed, in this new approach, it is possible that the 
underlying mechanism responsible for chaotic behavior can be characterized by only a 
few degrees of freedom. This new theory contradicts the earlier statistical theory 
which treats turbulence as a broad (i.e., many degrees of freedom) spectrum of 
periodic disturbances. But it is conjectured that it should include the secondary 
instability theory of transition. Perhaps, time-dependent instabilities are the 
likely mechanism by which an attractor remains strange. The arguments leading 
towards the concept of strange attractors and the Ruelle-Takens theory of turbulence 
will become evident in the section on results and their discussion. 

K. Ghia, Osswald and U. Ghia (1985a) analyzed massively separated flow past a 

12 percent thick Joukowski airfoil using a symmetric C-grid. The Reynolds number 
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ranged from 10 to 10 and the angle of attack a was varied up to 10°. For the 
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configuration of Re = 10 , aj.=15°, the unsteady massively separated flow asymptoted 
to a limit cycle. Subsequently, Osswald, K. Ghia and U. Ghia (1985b) significantly 
improved the above analysis by introducing circulation in the computation of 
conformal clustered C-grid. This resulted in the branch cut approaching the wake 
centerline a few chords downstream of the trailing edge (TE). A numerical method 
was developed to treat the branch cut of the C-grid implicitly while also 
appropriately treating the vorticity singularity at the TE. The versatility and 
efficiency of the numerical method were demonstrated by providing the massively 
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separated flow structure for Re=10 and a f =30° and also the transient results, at 
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early times, for Re=10 and a^=53°. The present study is a sequel to these earlier 
studies and has the following objectives: 

i. To provide the detailed flow structure of massively separated flow for 
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Re=10 , ot = 53°» in the post-stall regime of the flow past a Joukowski 
airfoil, 

ii. to obtain time-dependent aerodynamic lift, drag and moment coefficients 
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for Re=10 and between 15° to 53° and, finally, 

iii. to understand the observed quasiperiodicity and bring forth any similarity 
possible with strange attractors. 



ON THE ANALYSIS AND NUMERICAL METHODS 


The conservative form of the unsteady Navier-Stokes equations in terms of 
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vorticity and stream function in generalized curvilinear coordinates (£ , £ ) is 
used. The far-field boundary conditions of uniform flow are, strictly, valid only 
at infinity. To circumvent the large values of occurring in the far-field and 
also facilitate the numerical implementation of the far-field boundary conditions, 
the stream function i|> is decomposed into two parts such that = ^ in + ^ v > neither of 
which need to be small in comparison to the other. Here, is a solution of the 
stream function equation with zero vorticity and satisfies the free-stream condition 
at infinity and zero normal velocity at the surface of the airfoil. This 
decomposition is used so as to enable the far-field boundary condition to be placed 
at infinity for the viscous-flow calculations. 

The Reynolds number used in this study is defined as Re = U^c/v where c is the 
airfoil chord. The characteristic time is given as t = tVCc/U^), with being the 
undisturbed free-stream velocity. The boundary conditions correspond to uniform 
flow at infinity, together with the no-slip conditions along the body surface. The 
initial conditions for the viscous flow are taken as the corresponding steady-state 
inviscid solution. 

A clustered conformal grid is generated; the clustering is controlled by 

appropriate one-dimensional (1-D) stretching transformations. An attempt is made to 

resolve many of the multiple scales of the unsteady flow with massive separation, 

while maintaining the transformation metrics to be smooth and continuous in the 

entire flow field. Typical clustered conformal grids consisting of (230, 46) points 

for a 12 percent thick symmetric Joukowski airfoil with a = 15°, 30° and 53° are 
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shown in Fig. 1. Here, a is the grid angle of attack, i.e. the asymptotic slope of 
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the coordinate line emanating from the TE; this may, or may not, be the same as the 
flow angle of attack a f . 

A fully implicit time-marching method is developed such that all spatial 

derivatives are approximated using central differences and no use is made of any 

artificial dissipation. The numerical method solves the discretized equations using 

the alternatingj-direction implicit-block Gaussian elimination (ADI-BGE) method and 
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has overall 0[At, (A£ ) , (A £ ) ] accuracy. The boundary conditions on the 
vorticity are so implemented as to maintain second-order spatial accuracy. The 
results obtained for this study are discussed next. 

RESULTS AND DISCUSSION 

A 12 percent thick symmetric Joukowski airfoil is used in this study and has 
two especially attractive features, (i) The Joukowski airfoil can be accurately 
represented using conformal transformations; the details of these and the clustering 
transformations used were given by Osswald, K. Ghia and U. Ghia (1985a). (ii) The 
presence of a sharp TE leads to a much stronger interacting region and, hence, truly 


tests the analysis developed. This unsteady Naviei — Stokes analysis and the 

corresponding numerical method are used to study three flow configurations in 

detail. All of these configurations have the same Re = 1,000 but the value of the 

free-stream incidence varies such that a^. = 15°, 30° and 53 °, respectively. As 

stated earlier, a C-grid with (230x^6) mesh points has been used; this number was 

arrived at based on the available full-core capacity of the supermini Perkin Elmer 

3250 MPS computer system; this then also precluded a mesh refinement study. Hence, 

an attempt* was made to convert the program to run on the CYBER 205 VPS/32 computer 

at NASA Langley Research Center; this effort is still in progress. The CPU time t 

require to advance the solution by one time step per spatial grid point is referred 
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to as a "computational effort" index; the value achieved for t was 2.14x10 seconds 
for a Perkin Elmer 3250 MPS supermini computer system. 

Flow In Post-Stall Regime; Re » 1000, = 53° 

There are two critical parameters in the flow past an airfoil, namely, Re and 
a f . In the present study, Re is held constant and only is increased from 15° to 

53°, so that otj. becomes the critical parameter in this discussion. For the high 
angle-of-attack case with a^= 53°, the Joukowski airfoil appears, to the oncoming 
stream, as an apparent bluff body. The massively separated vortex-dominated flow in 
the post-stall regime for this configuration of the Joukowski airfoil, is 
exceedingly complex and, from the results available so far, it is feasible to 
conjecture two hypotheses; see Fig. 5. One possibility is that the solution has 
still not asymptoted to an exact limit cycle but may do so subsequently. The second 
possibility, based on the results for a f = 30 to be discussed later, is that the 
solution may asymptote to a quasiperiodic state, with anywhere from 3 to 8 
incommensurate frequencies. The second hypothesis is more likely to prevail and 
some arguments in its favor will be presented later. However, the discussion for 
this configuration is presently considered speculative in nature, since the results 
in Fig. 5 do not permit identifying a limit-cycle and describing its evolution. 
Hence, the time instants at which the results are presented, are selected 
arbitrarily to show the spatial structure between the LE-TE-LE sheddings. 

The instantaneous stream function contours presented in Figs. 2a, c, e, g show 
massively separated flow with large eddies present over the suction surface as well 
as in the wake. The stagnation streamline on the pressure surface fluctuates and is 
closer to LE when the LE eddy is stronger in strength (see Fig. 2a) as compared to 
Figs. 2c, e where the clockwise rotating eddy has already been shed and is in the 
process of strong interaction with the counterclockwise eddy from the TE. The 
corresponding vorticity contours are shown in Figs 2b, d, f, h and are of greater 


* This effort was made in August 1985, while the present investigators were in 
residence at ICASE. Subsequently, the computer program has been made fully 
operational and, it is anticipated thr.t the much-needed grid refinement in the 
normal direction will be carried out shortly. 



significance, since they are independent of the observer's reference frame and, 
hence, permit more meaningful assessment of the convection and diffusion of the 
vortices. Figure 2b shows the intense TE vortex and also the weaker LE vortex. As 
compared to the cases with a^. = 15° and 30°, the interaction between these two 
vortices is initially of a weaker nature. 

Figure 2d corresponds to intense interaction between the LE vortex and a part of the 
TE vortex, while the remaining part of the TE vortex is elongated and its shedding 
is imminent. Still there is a substantial strength elongated vortex that has 
remained as seen in Fig. 2f and it is about to shed at t = 52.7; also, the shear 
layer from the LE has thickened. Finally, in Fig. 2h, both the LE and TE vortices 
are beginning to intensify and the earlier interacting vortices are on the verge of 
separating from the LE shear layer. The subsequent flow structure at t = 54.5, not 
shown here, is very similar to that at t = 50 . -4 and corresponds to shedding of the 
vortex at the LE. This evolution process could have been visualized more clearly, 
had color contour plots been available to distinguish the LE and TE vortices. 

Limit Cycle Analysis 

i . Coefficients of Lift, Drag and Moment for Re ° 1000, a f = 15°, 30° and 53° 

Originally, the aerodynamic coefficients were computed only at intervals of 0.1 
characteristic time. This was quite satisfactory for qualitative assessment of the 
flow evolution but too coarse for its detailed analysis. Subsequently, the computer 
program has been modified to provide the coefficients of lift C^, drag and moment 
(nose-up positive about the quarter chord point) at every At increment. The 
configuration with = 15°, which requires minimum time to reach the time- 
asymptotic limit cycle solution has been completely recomputed, using a slightly 
improved grid near. the TE. The flow configuration with = 30° has been recomputed 
between t = 45 and t = 58, whereas, due to limitation of availability of CPU time on 
host computer system, the configuration with a^, = 53° is currently available with 
C , C and C computations at the interval of 0.1 t only. Figures 3a, b, c show C , 
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C D and C M corresponding to a f =15°. As seen in this figure, C^ rises initially but 
drops very sharply during the transient phase and asymptotes to a near-limit-cycle 
solution corresponding to the dominant frequency for the shedding of vortices from 
the TE. The L^ norm of the entire vorticity and stream function fields were 
carefully examined to ensure that the near-limit-cycle has been achieved. This 
limit-cycle solution is an "ordinary attractor", the attractor being the 1-D object 
to which the phase-space trajectories are attracted at all times. The complete 
motion is known once the geometry of the attractor is determined. Hence, it is also 
possible to compute the mean flow by averaging the flow over one complete cycle; 
similarly, it is possible to determine the Reynolds stresses from first principles, 
although these computations have not been performed in the present study. 

For the flow configuration with a f =30 o , the physics changes dramatically, as 
seen in Figs. 4a, b, c. The curves for the force coefficients show that one limit 


cycle consists of two TE vortices sheddings. There are now two shedding 
frequencies, or modes, associated with this more complex attractor. As shown in 
Fig. 4a, the first frequency is associated with the shedding which takes place at 
point 1, whereas the second frequency is associated with the shedding which takes 
place at point 3. The LE shear layer associated with the first frequency is thinner 
and more intense, as compared to that associated with the second frequency. The 
energy now oscillates between the two unstable modes through a nonlinear coupling. 
The appearance of subharmonics signifies small modulations in the shedding 
frequency. This flow field, with its two natural incommensurate frequencies, is 
referred to as a quas iper iodic flow, also known as "Hopf bifurcation" into an 
invariant torus. From Fig. 4a, it is clear that the initial state at point 10 of a 
new cycle is slightly different from that at point 8. If the phase-space 
trajectories were drawn, this solution may very well show the tendency to fill a 
rather significant surface area of the torus. Finally, it should be noted that the 
C D peaks in Fig. 4d correspond to points 2 and 9 in Fig. 4a. 

For the case with cij. = 53 0 » the results obtained up to t = 74 may be far from 
approaching an asymptotic state. Figure 5a shows some resemblance of quasiper iodic 
flow with three incommensurate frequencies, this fact being further supported by the 
curves in Figs. 5b, c. From their numerical experiments, Grebogi, Ott and 
Yorke (1 983) have also shown the existence of quas i per i od ic i ty with three 
incommensurate frequencies. The state of the system at a given time instant in one 
cycle is not quite repeated at the corresponding time instant in the subsequent 
cycle. The phase-space portrait, not shown here, is very complex, where the surface 
has folded repeatedly onto itself, so that it appears to be a strange attractor. 
This is an indication, although preliminary, that the flow may be exhibiting a route 
to chaos. There are a few rigorous approaches to characterize a strange attractor. 
These are the determination of (i) the Lyapunov exponent, (ii) the fractal dimension 
of the attractor, which is related to the number of degrees of freedom; and finally, 
(iii) the Kolmogorov entropy. These indices still need to be studied thoroughly in 
order to rigorously analyze the route to chaos in a meaningful way. 

i i . Entire State of the System 

Analysis of lift, drag and moment histories by themselves is not sufficient to 
examine limit-cycle behavior since it is quite possible that two different dynamic 
states of the system may produce the same values of C^, C D and C M « To avoid such a 
possibility, the norms of the difference between the evolving stream function and 
vorticity fields and their corresponding values at some designated initial time were 
computed at each time step. Only when this measure of the internal dynamic state of 
the system reached a value below a specified tolerance, was it considered that an 
initial state of the system had actually re-occurred. To date, this analysis has 
been carried out only for the configuration with Re=1 ,000 and a f =15°. Here, a near- 
limit-cycle (within the specified tolerance) does occur, with a period of 1.416 



characteristic-time units and a nondimensional shedding frequency of 0.706. The 
corresponding Strouhal number, S = fc sina/U^ = 0.18, (see K. Ghia, Osswald and 
U. Ghia (1985b)), agrees well with the universal wake-based number of Roshko (195^) . 
From Figs. 6a and c, as well as 6b and d, it is clear that the large as well as 
small-scale motion repeats itself and that a limit cycle has been nearly 
established. For the case of Re=1000 and 0^=30°, the L 2 -norm of the solution has 
yet to be examined. Instead, the limit cycle was arrived at from visual observation 
of the variation of aerodynamic coefficients presented in Fig. 4. As seen in Fig. 
4a, the two Strouhal numbers are S=0.15 and S=0.20, corresponding to two shedding 
frequencies. Considering the two sheddings per cycle, the overall Strouhal number S 
is obtained as 0.17 which, again, is in the correct range. The stream-function and 
vorticity contours in Figs. 7a-d for t=45.1 closely resemble those at 50.98, but the 
two states are not identical, implying a slightly different initial state each time. 
Finally, for Re=1000 and 0^=53°. the results obtained so far are not sufficient to 
warrant a limit-cycle analysis. It is likely that the flow field may asymptote to 
quasiperiodicity with three incommensurate frequencies, as can be inferred from Fig. 
5b. However, this is a very weak justification and, hence, no attempt is made to 
compute the Strouhal number yet. In order to understand the degree to which this 
motion repeats itself, points 1 and 14, corresponding to LE shedding after three 
sheddings, have been plotted in Figs. 8a-d at t=46.2 and 58.6. The stream-function 
and vorticity contours are similar, but there are sufficient differences. This 
leads to the conjecture that the motion may be quasiperiodic with three frequencies, 
and may be on the route to chaos. 

CONCLUSION 

The unsteady Navier-Stokes analysis is shown to be capable of analyzing the 
massively separated, persistently unsteady flow in the post-stall regime of a 
Joukowski airfoil for angle of attack as high as 53°. The analysis has provided the 
detailed flow structure, showing the complex vortex interaction for this 
configuration. The aerodynamic coefficients, C^» C D and C^ for lift, drag and 
moment, respectively, have been calculated; these aid further in characterizing this 
unsteady flow. The use of a f as the critical parameter appears to be appropriate, 
but further computations are needed to convincingly demonstrate the state of the 
flow for a f =53°. The phase-space portrait not shown for this case does, at least 
qualitatively show that the asymptotic solution is a strange attractor. This study 
has so far only computed the spatial structure of the vortex interaction. It is now 
important to potentially use ; the large-scale vortex interactions, an additional 
energy source, to improve the aerodynamic performance. To achieve this, a 
quantitative analysis must be carried out from the results obtained so far to 
analyze the dependence of vortex evolutions, their growth and strength, and the 
various interactions amongst vortices and with the airfoil. This latter study would 



aid in carefully characterizing the flow field as a function of the airfoil 
geometry. It is only after this phase that the undesirable aerodynamic features can 
be appropriately controlled and the energy of the unsteady vortex-dominated flow 
harnessed to improve the aerodynamic performance as desired under critical operating 
conditions. 
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Fig. 1. Typical Grid Distributions for Symmetric 
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Fig. 2. LE-TE-LE Vortex-Shedding Cycle for Joukowski Airfoil at Re = 1,000, 
ct^ = 53° . 











Fig. 4. Instantaneous Aerodynamic Coefficients at Re = 1,000, = 30°. 

(a) Lift C L , (b) Drag Cp, (c) Moment C^. 
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Fig. 8. LE-TE-LE Vortex Shedding for Joukowski Airfoil at Re = 1,000, = 53°. 
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